In this script, we created GAMs and GLMs to find out whether survival per stage or reproduction probabilities have a significant effect on lambda.

Setup

## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")

## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore 
## at a later stage
library("effects")
## Loading required package: carData
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library("gam")
## Loading required package: splines
## Loading required package: foreach
## Loaded gam 1.20
library("ggeffects")
library("here")
## here() starts at C:/Users/sinah/Documents/IZW/Drenske_2020_PVA_NBI
library("lattice")
library("mgcv")
## Loading required package: nlme
## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.
## 
## Attaching package: 'mgcv'
## The following objects are masked from 'package:gam':
## 
##     gam, gam.control, gam.fit, s
library("modelsummary")
library('MuMIn')
library("plot3D")
## Warning: package 'plot3D' was built under R version 4.0.5
library("pwr")
library("tidyverse")
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.0
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::collapse()   masks nlme::collapse()
## x dplyr::filter()     masks stats::filter()
## x dplyr::lag()        masks stats::lag()
## x purrr::when()       masks foreach::when()

Data

df_sc_mean <- readRDS(file = here("output", "data-proc", "model-proc", "df_sc_mean.rds"))

Management improvement scenarios

unique(df_sc_mean$Repro)
## [1] 0.53 0.58 0.66 1.06 1.41 3.97

Change the column names

colnames(df_sc_mean)[20:24] <- c("m1", "m2", "m3", "m4", "RR") 

Add columns for survival

df_sc_mean$s1 <- (1 - df_sc_mean$m1)
df_sc_mean$s2 <- (1 - df_sc_mean$m2)
df_sc_mean$s3 <- (1 - df_sc_mean$m3)
df_sc_mean$s4 <- (1 - df_sc_mean$m4)

Mortality and respective survival

print("m1 vs s1")
## [1] "m1 vs s1"
table(df_sc_mean$m1)
## 
##  0.2  0.3 0.36 
##  108  108  110
table(df_sc_mean$s1)
## 
## 0.64  0.7  0.8 
##  110  108  108
print("m2 vs s2")
## [1] "m2 vs s2"
table(df_sc_mean$m2)
## 
## 0.08 0.19 0.26 
##  108  108  110
table(df_sc_mean$s2)
## 
## 0.74 0.81 0.92 
##  110  108  108
print("m3 vs s3")
## [1] "m3 vs s3"
table(df_sc_mean$m3)
## 
## 0.14 0.24 0.31 
##  108  108  110
table(df_sc_mean$s3)
## 
## 0.69 0.76 0.86 
##  110  108  108
print("m4 vs s4")
## [1] "m4 vs s4"
table(df_sc_mean$m4)
## 
## 0.02 0.14 0.22 
##  108  108  110
table(df_sc_mean$s4)
## 
## 0.78 0.86 0.98 
##  110  108  108

Calculate extinction probability

df_sc_mean$pext <- df_sc_mean$Extinct/100

Explore database

names(df_sc_mean)
##  [1] "Scenario"        "Run"             "Year"            "IndAll"          "IndB"            "IndK"            "IndUb"           "IndFP"           "IndP"           
## [10] "IndBFP"          "IndBP"           "IndKFP"          "IndKP"           "IndUbFP"         "IndUbP"          "IndJuv"          "IndSub1"         "IndSub2"        
## [19] "IndAdu"          "m1"              "m2"              "m3"              "m4"              "RR"              "N_Supplements"   "Supplement_Time" "Stoch_Frequency"
## [28] "Stoch_Severity"  "N_Juveniles"     "N_Subadults1"    "N_Subadults2"    "N_Adults"        "Lambda_sto"      "SD_Lambda"       "Extinct"         "s1"             
## [37] "s2"              "s3"              "s4"              "pext"
View(df_sc_mean) # okay, these are mortalities.

unique(df_sc_mean[,c(20:24)]) # m1, m2, m3, m4, RR values per scenario
## # A tibble: 326 x 5
##       m1    m2    m3    m4    RR
##    <dbl> <dbl> <dbl> <dbl> <dbl>
##  1   0.2  0.08  0.14  0.02 0.53 
##  2   0.2  0.08  0.14  0.02 0.580
##  3   0.2  0.08  0.14  0.02 0.66 
##  4   0.2  0.08  0.14  0.02 1.06 
##  5   0.2  0.08  0.14  0.14 0.53 
##  6   0.2  0.08  0.14  0.14 0.580
##  7   0.2  0.08  0.14  0.14 0.66 
##  8   0.2  0.08  0.14  0.14 1.06 
##  9   0.2  0.08  0.14  0.22 0.53 
## 10   0.2  0.08  0.14  0.22 0.580
## # ... with 316 more rows
unique(df_sc_mean$s1)
## [1] 0.80 0.70 0.64
unique(df_sc_mean$s2)
## [1] 0.92 0.81 0.74
unique(df_sc_mean$s3)
## [1] 0.86 0.76 0.69
unique(df_sc_mean$s4)
## [1] 0.98 0.86 0.78
unique(df_sc_mean$RR)
## [1] 0.53 0.58 0.66 1.06 1.41 3.97
# start with pExt
table(df_sc_mean$pext)#df_sc_StevS_mean
## 
##    0 0.01 0.02 0.03 0.06 0.24 
##  316    4    3    1    1    1
#  0    0.01 0.02 0.03 0.06 0.24 
# 316    4    3    1    1    1 

# Comment: Only in 10 out of 326 scenarios did the population ever go extinct. 
# With such a distribution of the response variable, 
# there is no need to calculate a model. 
# Therefore we have described the distribution of survival and reproduction values

rows_ext_scen <- which(df_sc_mean$pext > 0)
ext_scen <- df_sc_mean[rows_ext_scen,]
ext_scen
## # A tibble: 10 x 40
##    Scenario    Run  Year IndAll  IndB  IndK IndUb IndFP  IndP IndBFP IndBP IndKFP IndKP IndUbFP IndUbP IndJuv IndSub1 IndSub2 IndAdu    m1    m2    m3    m4    RR N_Supplements
##    <fct>     <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>   <dbl>  <dbl>  <dbl>   <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>         <dbl>
##  1 0.20.26~ 10451.  24.9   47.4     0     0     0     0  41.3      0     0      0     0       0      0  0.728   11.7     8.67   26.3  0.2   0.26  0.31  0.22 0.53              0
##  2 0.30.19~ 17651.  25.0   42.7     0     0     0     0  36.6      0     0      0     0       0      0  0.727    9.65    7.83   24.5  0.3   0.19  0.31  0.22 0.53              0
##  3 0.30.26~ 21250.  24.9   31.7     0     0     0     0  25.8      0     0      0     0       0      0  0.729    7.37    5.51   18.1  0.3   0.26  0.31  0.22 0.53              0
##  4 0.30.26~ 21351.  25.0   41.4     0     0     0     0  35.5      0     0      0     0       0      0  0.727    9.98    7.43   23.3  0.3   0.26  0.31  0.22 0.580             0
##  5 0.360.1~ 27250.  25     42.9     0     0     0     0  36.7      0     0      0     0       0      0  0.725    9.18    7.46   25.5  0.36  0.19  0.24  0.22 0.53              0
##  6 0.360.1~ 28451.  24.9   31.9     0     0     0     0  26.0      0     0      0     0       0      0  0.730    6.90    5.71   18.6  0.36  0.19  0.31  0.22 0.53              0
##  7 0.360.1~ 28551.  24.9   42.1     0     0     0     0  36.1      0     0      0     0       0      0  0.727    9.52    7.76   24.1  0.36  0.19  0.31  0.22 0.580             0
##  8 0.360.2~ 30851.  24.9   34.2     0     0     0     0  28.2      0     0      0     0       0      0  0.729    7.51    5.69   20.2  0.36  0.26  0.24  0.22 0.53              0
##  9 0.360.2~ 50051.  24.3   21.4     0     0     0     0  15.6      0     0      0     0       0      0  0.750    4.84    3.54   12.3  0.36  0.26  0.31  0.22 0.53              0
## 10 0.360.2~ 32151.  25.0   29.8     0     0     0     0  24.1      0     0      0     0       0      0  0.727    6.90    5.21   17.0  0.36  0.26  0.31  0.22 0.580             0
## # ... with 15 more variables: Supplement_Time <dbl>, Stoch_Frequency <dbl>, Stoch_Severity <dbl>, N_Juveniles <dbl>, N_Subadults1 <dbl>, N_Subadults2 <dbl>, N_Adults <dbl>,
## #   Lambda_sto <dbl>, SD_Lambda <dbl>, Extinct <int>, s1 <dbl>, s2 <dbl>, s3 <dbl>, s4 <dbl>, pext <dbl>
table(ext_scen$m1)
## 
##  0.2  0.3 0.36 
##    1    3    6
table(ext_scen$m2)
## 
## 0.19 0.26 
##    4    6
table(ext_scen$m3)
## 
## 0.24 0.31 
##    2    8
table(ext_scen$m4)
## 
## 0.22 
##   10
table(ext_scen$RR)
## 
## 0.53 0.58 
##    7    3
# We can see that the highest mortality of adults was part of all scenarios where the population went extinct
# See also the excel file under "output\data-proc\model-proc\08_PVA_GLM_table_ext_scen.xlsx"
# Explore lambda for s1 and RR
x <- df_sc_mean$s1
y <- df_sc_mean$RR
z <- df_sc_mean$Lambda_sto
scatter3D(x, y, z, pch = 19, cex = 0.7, clab = c("Lambda"), xlab = "s1", ylab = "RR", zlab = "Lambda")

contourplot(Lambda_sto ~ s1 * RR, data = df_sc_mean, labels = FALSE,
            region = TRUE)

plot(df_sc_mean$Lambda_sto)
abline(h=1) # Lambda > 1 indicates population growth

# Explore lambda for s4 and RR
x <- df_sc_mean$s4
y <- df_sc_mean$RR
z <- df_sc_mean$Lambda_sto
scatter3D(x, y, z, pch = 19, cex = 0.7, clab = c("Lambda"), xlab = "s4", ylab = "RR", zlab = "Lambda")

contourplot(Lambda_sto ~ s4 * RR, data = df_sc_mean, labels = FALSE,
            region = TRUE)

GAM

A general additive model (GAM) tests for non-linearity

The extinction probability is the response variable and the predictor variables are formed by the survival probabilities per stage (s1-s4) and the reproduction (RR)

# define a k
# myk - stands for 'my number of knots/ smoothing splines' :-)
myk <- 3 # or 4

Response = Lambda

## Lambda and survival
gammod_2_1 <- gam(formula = Lambda_sto ~ s(s1, k=myk) + s(s2, k=myk) + 
                          s(s3, k=myk) + s(s4, k=myk) +  
                          s(RR, k=myk), 
                          data = df_sc_mean, 
                          na.action = 'na.fail', family=Gamma) 
summary(gammod_2_1) 
## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Lambda_sto ~ s(s1, k = myk) + s(s2, k = myk) + s(s3, k = myk) + 
##     s(s4, k = myk) + s(RR, k = myk)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9036404  0.0004053    2229   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df      F p-value    
## s(s1) 1.000  1.000  964.7  <2e-16 ***
## s(s2) 1.000  1.000  887.0  <2e-16 ***
## s(s3) 1.000  1.000  893.8  <2e-16 ***
## s(s4) 1.908  1.991 5617.9  <2e-16 ***
## s(RR) 1.999  2.000 5328.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.985   Deviance explained = 98.6%
## GCV = 6.6955e-05  Scale est. = 6.5312e-05  n = 326
# F values for s4 and RR are very high - they have the highest influence

## Lambda and mortality
gammod_2_2 <- gam(formula = Lambda_sto ~ s(m1, k=myk) + s(m2, k=myk) + 
                          s(m3, k=myk) + s(m4, k=myk) +  
                          s(RR, k=myk), 
                          data = df_sc_mean, 
                          na.action = 'na.fail', family=Gamma) 
summary(gammod_2_2) 
## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Lambda_sto ~ s(m1, k = myk) + s(m2, k = myk) + s(m3, k = myk) + 
##     s(m4, k = myk) + s(RR, k = myk)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9036404  0.0004053    2229   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df      F p-value    
## s(m1) 1.000  1.000  964.7  <2e-16 ***
## s(m2) 1.000  1.000  887.0  <2e-16 ***
## s(m3) 1.000  1.000  893.8  <2e-16 ***
## s(m4) 1.908  1.991 5617.9  <2e-16 ***
## s(RR) 1.999  2.000 5328.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.985   Deviance explained = 98.6%
## GCV = 6.6955e-05  Scale est. = 6.5312e-05  n = 326
# F values for s4 and RR are very high - they have the highest influence

# no difference between survival and mortality

Plots

plot(gammod_2_1)

vis.gam(gammod_2_1, view = c("s4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
        ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
        xlab = "s4", ylab = "RR", zlab = "Lambda", 
        color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)

vis.gam(gammod_2_2, view = c("m4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
        ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
        xlab = "m4", ylab = "RR", zlab = "Lambda", 
        color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)

Save the plot

# # s4 and RR
# png(filename = here("plots", "04_PVA", "persp_s4_RR_lambda_20211124.png"), width = 1200, height = 772)
# vis.gam(gammod_2_1, view = c("s4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
#         ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
#         xlab = "s4", ylab = "RR", zlab = "Lambda", 
#         color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)
# dev.off()

Effect size

Calculating the effect size

pwr.t.test(n=326, sig.level=0.05, power = 0.8, alternative="greater")
## 
##      Two-sample t test power calculation 
## 
##               n = 326
##               d = 0.1949595
##       sig.level = 0.05
##           power = 0.8
##     alternative = greater
## 
## NOTE: n is number in *each* group

GLM

Response = Lambda

# With interactions between RR and all survival values
mod1_3_1 <- glm(Lambda_sto ~ (s1 + s2 + s3 + s4) + RR, 
              data = df_sc_mean, na.action = 'na.fail', family=Gamma)

# With interactions between RR and all survival values
mod1_3_2 <- glm(Lambda_sto ~ (m1 + m2 + m3 + m4) + RR, 
              data = df_sc_mean, na.action = 'na.fail', family=Gamma)
summary(mod1_3_1)  
## 
## Call:
## glm(formula = Lambda_sto ~ (s1 + s2 + s3 + s4) + RR, family = Gamma, 
##     data = df_sc_mean, na.action = "na.fail")
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.174754  -0.010237  -0.002538   0.005698   0.039817  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.870086   0.018618  100.45   <2e-16 ***
## s1          -0.199206   0.012780  -15.59   <2e-16 ***
## s2          -0.170727   0.011390  -14.99   <2e-16 ***
## s3          -0.182280   0.012103  -15.06   <2e-16 ***
## s4          -0.522307   0.010212  -51.15   <2e-16 ***
## RR          -0.122401   0.002476  -49.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.000289145)
## 
##     Null deviance: 1.507876  on 325  degrees of freedom
## Residual deviance: 0.095412  on 320  degrees of freedom
## AIC: -1645.8
## 
## Number of Fisher Scoring iterations: 4
anova(mod1_3_1, test='LRT') 
## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Lambda_sto
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   325    1.50788              
## s1    1  0.05411       324    1.45377 < 2.2e-16 ***
## s2    1  0.04857       323    1.40520 < 2.2e-16 ***
## s3    1  0.04808       322    1.35712 < 2.2e-16 ***
## s4    1  0.67890       321    0.67822 < 2.2e-16 ***
## RR    1  0.58281       320    0.09541 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# s4 and RR have the highest deviance - they explain the model the most

#dredge(mod1_3_1)
plot(allEffects(mod1_3_1),type='response', ylim = c(1,2))

Save the plot

png(filename = here("plots", "04_PVA", "lambda_effects_plot20211216.png"), width = 1200, height = 772)
plot(allEffects(mod1_3_1),type='response', ylim = c(1,2))
dev.off()
## png 
##   2

R square R^2= 1-(Residual deviance/Null deviance)

round(1 - (0.095412/1.507876), digits = 2) 
## [1] 0.94
summary(mod1_3_2) 
## 
## Call:
## glm(formula = Lambda_sto ~ (m1 + m2 + m3 + m4) + RR, family = Gamma, 
##     data = df_sc_mean, na.action = "na.fail")
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.174754  -0.010237  -0.002538   0.005698   0.039817  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.795567   0.005226  152.23   <2e-16 ***
## m1           0.199206   0.012780   15.59   <2e-16 ***
## m2           0.170727   0.011390   14.99   <2e-16 ***
## m3           0.182280   0.012103   15.06   <2e-16 ***
## m4           0.522307   0.010212   51.15   <2e-16 ***
## RR          -0.122401   0.002476  -49.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.000289145)
## 
##     Null deviance: 1.507876  on 325  degrees of freedom
## Residual deviance: 0.095412  on 320  degrees of freedom
## AIC: -1645.8
## 
## Number of Fisher Scoring iterations: 4
anova(mod1_3_2, test='LRT') 
## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Lambda_sto
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   325    1.50788              
## m1    1  0.05411       324    1.45377 < 2.2e-16 ***
## m2    1  0.04857       323    1.40520 < 2.2e-16 ***
## m3    1  0.04808       322    1.35712 < 2.2e-16 ***
## m4    1  0.67890       321    0.67822 < 2.2e-16 ***
## RR    1  0.58281       320    0.09541 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# s4 and RR have the highest deviance - they explain the model the most

#dredge(mod1_3_2)
plot(allEffects(mod1_3_2),type='response')

R square R^2= 1-(Residual deviance/Null deviance)

round(1 - (0.095412/1.507876), digits = 2) 
## [1] 0.94

Model summary

datasummary_skim(df_sc_mean)
Unique (#) Missing (%) Mean SD Min Median Max
Run 326 0 16525.0 10085.5 50.5 16300.5 60150.5
Year 207 0 19.5 5.6 6.6 20.3 25.0
IndAll 326 0 830.9 473.9 21.4 1162.3 1406.3
IndB 1 0 0.0 0.0 0.0 0.0 0.0
IndK 1 0 0.0 0.0 0.0 0.0 0.0
IndUb 1 0 0.0 0.0 0.0 0.0 0.0
IndFP 1 0 0.0 0.0 0.0 0.0 0.0
IndP 326 0 813.2 465.2 15.6 1143.5 1386.3
IndBFP 1 0 0.0 0.0 0.0 0.0 0.0
IndBP 1 0 0.0 0.0 0.0 0.0 0.0
IndKFP 1 0 0.0 0.0 0.0 0.0 0.0
IndKP 1 0 0.0 0.0 0.0 0.0 0.0
IndUbFP 1 0 0.0 0.0 0.0 0.0 0.0
IndUbP 1 0 0.0 0.0 0.0 0.0 0.0
IndJuv 193 0 1.0 0.4 0.7 0.9 2.6
IndSub1 326 0 214.3 135.0 4.8 252.5 674.8
IndSub2 326 0 152.1 91.8 3.5 181.6 351.3
IndAdu 326 0 463.5 266.3 12.3 576.7 787.2
m1 3 0 0.3 0.1 0.2 0.3 0.4
m2 3 0 0.2 0.1 0.1 0.2 0.3
m3 3 0 0.2 0.1 0.1 0.2 0.3
m4 3 0 0.1 0.1 0.0 0.1 0.2
RR 6 0 0.7 0.3 0.5 0.7 4.0
N_Supplements 1 0 0.0 0.0 0.0 0.0 0.0
Supplement_Time 1 0 4.0 0.0 4.0 4.0 4.0
Stoch_Frequency 1 0 0.0 0.0 0.0 0.0 0.0
Stoch_Severity 1 0 0.2 0.0 0.2 0.2 0.2
N_Juveniles 1 0 37.0 0.0 37.0 37.0 37.0
N_Subadults1 1 0 11.0 0.0 11.0 11.0 11.0
N_Subadults2 1 0 8.0 0.0 8.0 8.0 8.0
N_Adults 1 0 18.0 0.0 18.0 18.0 18.0
Lambda_sto 326 0 1.1 0.1 1.0 1.1 1.4
SD_Lambda 326 0 0.0 0.0 0.0 0.0 0.0
Extinct 6 0 0.1 1.4 0 0.0 24
s1 3 0 0.7 0.1 0.6 0.7 0.8
s2 3 0 0.8 0.1 0.7 0.8 0.9
s3 3 0 0.8 0.1 0.7 0.8 0.9
s4 3 0 0.9 0.1 0.8 0.9 1.0
pext 6 0 0.0 0.0 0.0 0.0 0.2
varnam <- c('s1' = 'Survival of stage 1',
            's2' = 'Survival of stage 2',
            's3' = 'Survival of stage 3',
            's4' = 'Survival of stage 4',
            'RR' = 'Reproductive Rate')
varnam2 <- c('m1' = 'Mortality of stage 1',
            'm2' = 'Mortality of stage 2',
            'm3' = 'Mortality of stage 3',
            'm4' = 'Mortality of stage 4',
            'RR' = 'Reproductive Rate')


# modelplot(mod1_3,
#           coef_omit = 'Interc',
#           coef_map = varnam,
#           conf_level = .99) 

modelplot(mod1_3_1, 
          coef_omit = 'Interc',
          coef_map = varnam, 
          conf_level = .99) 

modelplot(mod1_3_2, 
          coef_omit = 'Interc',
          coef_map = varnam2, 
          conf_level = .99) 

# You can see that s4/m4 and RR are important and s1-s3/m1-m3 are not.

Session Info

## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()
## [1] "2022-01-06 18:55:07 CET"
#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4        readr_1.4.0        tidyr_1.1.2        tibble_3.0.3       ggplot2_3.3.2      tidyverse_1.3.0   
## [10] pwr_1.3-0          plot3D_1.4         MuMIn_1.43.17      modelsummary_0.6.5 mgcv_1.8-33        nlme_3.1-149       lattice_0.20-41    here_0.1           ggeffects_1.0.0   
## [19] gam_1.20           foreach_1.5.1      effects_4.2-0      carData_3.0-4     
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4       colorspace_1.4-1  ellipsis_0.3.1    showtext_0.9      sjlabelled_1.1.7  rprojroot_1.3-2   estimability_1.3  fs_1.5.0          rstudioapi_0.11  
## [10] farver_2.0.3      showtextdb_3.0    remotes_2.2.0     fansi_0.4.1       lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  knitr_1.30        pkgload_1.1.0    
## [19] jsonlite_1.7.1    nloptr_1.2.2.2    broom_0.7.4       dbplyr_1.4.4      compiler_4.0.3    d6_0.1.0.0        httr_1.4.2        backports_1.1.10  assertthat_0.2.1 
## [28] Matrix_1.2-18     survey_4.0        cli_2.0.2         htmltools_0.5.0   prettyunits_1.1.1 tools_4.0.3       misc3d_0.9-0      gtable_0.3.0      glue_1.4.2       
## [37] tables_0.9.6      tinytex_0.26      Rcpp_1.0.5        cellranger_1.1.0  vctrs_0.3.4       iterators_1.0.13  insight_0.11.1    xfun_0.18         ps_1.4.0         
## [46] testthat_2.3.2    lme4_1.1-25       rvest_0.3.6       lifecycle_0.2.0   devtools_2.3.2    statmod_1.4.35    MASS_7.3-53       scales_1.1.1      hms_0.5.3        
## [55] yaml_2.2.1        memoise_1.1.0     stringi_1.5.3     highr_0.8         desc_1.2.0        checkmate_2.0.0   boot_1.3-25       pkgbuild_1.1.0    rlang_0.4.7      
## [64] pkgconfig_2.0.3   evaluate_0.14     labeling_0.3      processx_3.4.4    tidyselect_1.1.0  magrittr_1.5      bookdown_0.20     R6_2.4.1          generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       withr_2.3.0       survival_3.2-7    nnet_7.3-14       modelr_0.1.8      crayon_1.3.4      utf8_1.1.4       
## [82] rmarkdown_2.4     sysfonts_0.8.1    usethis_1.6.3     grid_4.0.3        readxl_1.3.1      blob_1.2.1        callr_3.5.0       rmdformats_0.3.7  webshot_0.5.2    
## [91] reprex_0.3.0      digest_0.6.25     stats4_4.0.3      munsell_0.5.0     viridisLite_0.3.0 kableExtra_1.3.1  mitools_2.4       tcltk_4.0.3       sessioninfo_1.1.1